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FOREWORD 


Occasionally  the  engineer  or  analyst  is  confronted  with  the 
problem  of  dealing  with  a  collection  of  moving  particles.  Fragments 
from  exploding  bombs,  water  from  a  sprinkling  system,  clouds  of 
meteoroids  in  the  vicinity  of  the  earth,  photons  from  a  distant  star, 
cosmic  rays  interacting  with  the  magnetic  field  of  the  earth,  all  have 
a  common  thread  in  theoretical  mechanics  and  are  connected  together  in 
this  report.  Not  much  of  the  material  is  new,  but  the  collection, 
simplification  and  unification  of  work  from  several  specialized  fields 
seemed  worthwhile. 
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1.  INTRODUCTION 


Imagine  a  swarm  of  independent  particles  moving  in  space  with  the 
ith  particle  having  a  location  defined  by  the  position  vector  r^  and  a 

velocity  v^.  A  unit  volume  with  its  center  located  by  the  position 

vector  r  will  contain  N(r,t)  particles  at  time  t.  The  number  N(r,t), 
called  the  particle  density,  may  fluctuate  in  time  and  may  be  viewed  as 
a  probability  density  function  as  is  common  in  statistical  or  quantum 
mechanics.  The  expression. 


d  N(r,t)  =  N(r,t)  dt. 


1.1 


may  be  defined  as  the  number  of  particles  expected  in  an  element  of 
volume  dr.  In  an  analogous_manner,  a  swarm  of  particles  can  be  described 
by  a  density  function  N(r,p,t)_  in  the  six  dimensional  space  of  position 
vector  r  and  momentum  vector  p  such  that  the  expression. 


d  N(r,p,t)  =  N(r,p,t)  dx  dy  dz  dpxdpydpz  , 


1.2 


is  the  probable  number  of  particles  contained  at  time  t  in  the  element 


of  volume  dxdydzdp^dp^dp^ . 


By  definition  the  expression. 


N(r,t) 


N(r,p,t) 


dPxdPydPz> 


1.3 


defines  the  density  of  particles  in  the  three-space  defined  by  the 
coordinates  x,  y  and  z. 


If  the  individual  particles  are  moving  in  a  force  field  such  that  a 
Hamiltonian  H  can  be  written  for  ea£h_of  them,  a  Hamiltonian_for  the 
swarm  can  be  written  in  terms  of  N(r,p,t)  and  the  variables  r  and  p 
such  that 


and 


3H/3xi 


-Pj  =  dp./dt 


3H/3p^  =  jL  =  dx^/dt 


1.4 


1.5 


where  the  x^  are  the  position  coordinates  x,  y  and  z  and  the  p^  are 

the  momenta  associated  with_them.  An  arbitrary  volume  in  the  six 
dimensional  space  of  r  and  p  can  be  defined  by  particles  on  its  surface 
as 


V 


dxjdx2dx3dpjdp2dp.. 
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so  that  from  the  divergence  theorem. 


av/3t  -////If’ 


v  dx, 


1.7 


where  V  •  v  is  defined  in  the  six  space  of  r  and  p  by  the  Einstein 
summation  convention  as 


V  •  v  dt  =  [3)L/3x^  +  3p^/3p^]dx 


1.8 


Because  of  Eqs  1.4  and  1.5,  Eq  1.8  becomes 
V  •  v  d: 


'  32H 

3x.  3p. 
1 


32H 


3p.  3x. 

1  l  l 


dt  =  0  *dx  . 


Thus , 


1.9 


dV/dt  =  0 


for  the  volume  defined  in  6-space  by  a  surface  of  moving  particles. 
In  particular,  the  elemental  volume  dx_,  associated  with  a  particular 
particle  as  it  moves  in  the  space  of  r  and  v,  is  a  constant  of  the 
motion.  Thus,  a  collection  of  particles  contained  in  an  element  of 
volume  dx  may  be  defined  at  points  1  and  2  as 


N(r1,p1,t1)  dij  =  N(r2,p2,t2)dx2.  1.10 

Since  dx^  and  dx2  are  equal, 

NCrj.Pj.tj)  =  N(r2,p2,t2),  1.11 

and  N(r,p,t)  is  also  seen  to  be  a  constant  along  a  particle  trajectory 
in  the  space  of  r  and  p.  This  particular  form  of  tht  Liouville  theorem 
has  proved  useful  in  the  treatment  of  meteoroid  distributions  and 
collisionless  plasmas  interacting  with  the  magnetic  field  of  the  earth. 
These  particular  applications  will  be  discussed  more  in  detail  later. 

A  detector  will  be  defined  as  any  surface  sensitive  to  particle 
impacts.  For  a  planar  surface  with  a  normal  ri,  the  rate  of  impact 
<j>(r,ri)  from  one  side  is  given  by  the  integral. 


<f>(r,n) 


N(r,v) 


v  •  ri  dv  dS, 


1.12 
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where  N(r,v)  represents  the  distribution  of  number  with  respect  to 
position  vector  r  and  velocity  v,  dv  is  the  element  of  volume  in  velocity 
space,  and  dS  is  an  element  of  area.  If  the  area  is  small  and  equal  to 
A,  Eq  1.12  can  be  written  as 


N(r,v)  v  •  n  dv 


The  simplest  example  of  particle  flow  is  that  of  a  uniform  monoenergetic 
monodirectional  stream  of  particles  with_density  N  and  velocity  V. 

In  this  case,  Eq  1.13  becomes  cf)(r,n)  =  NV*nA  =  JA  cos9 


where 


J  =  Nv 


cose  =  v  •  n/v. 


For  a  uniform  monoenergetic  isotropic  distribution  of  speed  v, 

N(r,v)  =  N/4tt 
v  •  n  s  vcos9 

and  for  a  planar  detector  of  unit  area. 


4>(r,n)  =  (N/4tt) 


'  0=ir/2 


v cos 0  sine  de  dcf>  =  Nv/2  . 


It  is  common  in  particle  detection  to  define  an  isotropic  spherical 
detector  which  presents  unit  area  to  all  directions.  The  count  rate 
<t>  for  such  a  detector  is  given  by 


=  J  Jn(t, 


v)  v  dv  • 


For  the  uniform  monoenergetic  monodirectional  stream  above. 


and  for  the  isotropic  distribution, 

4>  =  Nv.  1.19 

The  measurement  <J>  described  by  Eq  1.17  is  usually  called  the  particle  flux 
and  is  measured  in  p/m^  sec.  In  a  unit  volume,  if  there  are  particles 

with  speed  V^,  <t>  can  be  defined  as 


J 


1.20 


where 


and 


<P  =  Z  N.v.  =  Nv 
.11 
1 

N  =  Z  N. 

1 

1 


v-Z  Nivi/Z  Ni . 

i  i 


1.21 


1.22 


Occasionally  we  need  to  compute  the  rate  j£(r,h)  at  which  particles 
fall  on  an  element  of  surface  n  ds  located  at  r  when  the  particle 
source  is  located  at  the  coordinate  origin  and  has  an  emission 
described  by  S(V)  p/ster  sec.  In  this  case 


Hr,n)  =  s(V)  ds  (r*n)/r2.  1.23 

2 

The  expression  ds  r-n/r  is  defined  as  the  solid  angle  subtended  by 
the  surface  ds  as  measured  from  the  source  located  at  the  origin. 


A  particularly  interesting  and  useful  interpretation  of  the 
Liouville  theorem  has  been  applied  to  the  motion  of  charged  particles 
in  electromagnetic  fields,  a  subject  of  interest  to  geophysicists 
studying  the  population  of  particles  captured  in  the  magnetic  field 
of  the  earth  and  to  particle  physicists  trying  to  promote  thermonuclear 
reactions  among  particles  contained  in  a  magnetic  bottle.  Suppose  that 
a  particle  trajectory  is  defined  by  r(t)  and  that  its  velocity  is 
therefore  7(r,t)  =  dr/dt.  Suppose  that,  at  some  point  r,  on  the 
particle  trajectory,  a  directional  counter  sensitive  to  an  element  of 
solid  angle  dw  is^  pointed  along  the  particle  trajectory  so  that  its 
counting  rate  C(r^)  is  given  by 


Cfr^Vj) 


NOj.v) 


vdvdw  =  dw 


N(r^,Vj)  vdv 


1.24 


Since  a  constant  magnetic  field  does  not  change  the  speed  of  a  particle, 
at  some  other  point  r^  on  the  particle  trajectory. 


c(r2>v2)  =  CCrj.Vj)  . 


1.25 


where  it  is  understood  that  both  counting  rates  are  obtained  by  pointing 
the  directional  detector  along  the  trajectory  of  a  particular  particle. 
Eq  1.25  is  just  an  experimental  way  of  expressing  the  Liouville  theorem 
for  a  limited  application.  The  integral  in  Eq  1.24  gives  the  counting 
rate  for  a  directional  detector  but  also  defines  a  small  element  of 
volume  in  r-p  space.  When  this  same  element  of  volume  in  phase  space  is 
sampled  at  another  point  (or  an  equivalent  volume  in  a  time  independent 
system  of  particles),  it  must  contain  the  same  number  of  particles  and 
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result  in  the  same  counting  rate. 

When  the  first  American  satellites  were  placed  in  orbit  about  the 
earth,  they  carried  directional  detectors  for  charged  particles.  It  was 
possible  by  means  of  Eq  1.25  to  infer  the  particle  density  for  the 
entire  region  around  the  earth  by  making  measurements  only  in  the 
equatorial  plane. 


2.  COORDINATE  TRANSFORMATIONS  FOR  BASIC  PHYSICAL  QUANTITIES 

Occasionally  it  is  easier  to  perform  an  analysis  in  one  coordinate 
system  and  to  transform  to  another  for  which  an  answer  is  needed  than 
to  perform  an  analysis  in  the  final  coordinate  system.  Also, 
measurements  may  be  made  in  the  wrong  coordinate  system  for  easy 
application.  If  a  particle  has  a  velocity  v  in  the  primary  reference 
system  and  a  velocity  v"  in  another  system  moving  at  a  velocity  V 
relative  to  the  primary  system,  the  velocities  are  related  by  the 
classical  addition  law, 


v  =  V  +  v',  2.1 

and  the  position  vectors  are  related  by  the  equation, 

_  r  =  R  +  r',  2.2 

where  R  connects  the  coordinate  origins. 

If  a  collection  of  particles  are  observed  in  each  system  and  the  total 
number  of  particles  N  in  an  arbitrary  volume  is  given  by  the  equations 


N  =  JJ I j J  j  N(r,v,t)  dx  N'(r',  v',t)  dx', 

the  relationship  of  the  density  functions  is  given  by 


2.3 


so  that 


N(r,v,t)  dx  =  N[r(r')v(v')]  D(*j* V— -)  dx 

r ',  v ' 


2.4 


N'Cr'.v")  =  N[r(r'l,  v(v')]  ).  2.5 

r',v' 

The  symbol  D  is  the  Jacobian  associated  with  the  transformation  of 
coordinate  differentials. 


As  an  example,  suppose  that  the  distribution  in  the  first  system_ 
is  isotropic  and  that  the  distribution  in  a  system  moving  with  speed  V 
along  the  zz'  avis  is  desired.  We  will  define 
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2 

dt  =  dx  dy  dz  V  sine  dvde  d<f> 


2.6 


in  the  first  system.  Note  that  cartesian  coordinates  are  used  in  r 
space  and  spherical  coordinates  are  used_in  v  space.  If  the  coordinate 
origins  and  coordinate  axes  coincide  in  r  space,  from  Eq  2.2 


and 

From  Eq  2.1 

and 


r  =  r' 

dx  dy  dz  =  dx'  dy'  dz'. 

v  =  V  +  7' 

v2  =  V2  +  (v')2  +  2Vv'  cose. 


2.7 

2.8 


2.9 

2.10 


If  the  coordinate  areas  are  chosen  so  that  V  coincides  with  the  z  and  z' 


axes,  e  becomes  0'  in  the  r'0'<j>'  system,  as  shown  in  Figure  1.  It  is 
evident  from  Figure  1  that 

v2  =  V2  +  Cv')2  +  2  W'cos0'  2.  II 

V'  sin6'  =  v  sin0  2. 12 

<T  =  <f>  2.13 


so  that  the  Jacobian  can  be  written  after  some  algebra  as 


D(I^_)  = 

r',v' 


Thus 

N'(r.'v')  =  N [r (r ') ,  v(v')]v'/v.  2.15 

In  the  limit  as  v'  becomes  small,  v  and  V  become  the  same,  as  seen 
from  Figure  1  or  Eq  2.11,  and  the  distribution  becomes  monoenergetic 
and  monodirectional  with  velocity  V.  If  V  becomes  small,  v  and  v' 
become  the  same  and  the  distributions  become  equal.  For  the  particular 
case  of  a  uniform  isotropic  distribution  such  that 

N  (r, v)  =  A  ,  Vj<v<v2,  2.16 

the  particles  are  contained  in  velocity  space  between  two  spheres  of 


3v 

30 

3v ' 

3v  ' 

3v 

30 

30' 

30' 
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radius  and  v2>  The  density,  4tt  A(v^  -  v^)/3,  is  just  the  volume 

between  these  two  spheres  multiplied  by  A.  In  the  primed  coordinate 
system,  the  use  of  Eq  2.10  shows  the  particles  to  be  contained  in 
the  volumes  between  the  surfaces  given  by  the  equations, 

V2  +  (v')2  +  2  Vv'  cos0'  =  v22 

and 

V2  +  O')2  +  2  Vv'  cos9'  =  Vj2. 

These  are  spheres  with  centers  located  at  -V  on  the  polar  axis  in 
velocity  space  and  with  radii  v2  and  v^.  Thus  the  particle  density 

is  conserved  by  the  transformation  and  the  distribution  is  changed 
according  to  Eq  2.15.  Eq  2.15  can  be  integrated  directly  over  velocity 
space  to  obtain  the  particle  density,  but  considerable  algebra  is 
involved.  As  V  becomes  large,  Eqs  2.17  and  2.18  show  that  the  velocity 
vectors  for  particles  lie  between  the  cones  defined  by  the  inequalities, 

it  -sin  1  Oj/V)  <_  0'  <_  tt  -sin  *(v2/V).  2.19 

In  the  limit  as  V  becomes  much  greater  than  v^  and  v2,  an  observer 
in  the  primed  system  sees  all  the  particles  traveling  with  velocity  -V. 

In  classical  mechanics,  there  are  a  number  of  quantities  which  are 
independent  of  the  frame  of  reference  and  are  said  to  be  invariant  to 
translational  and  rotational  coordinate  transformations.  These 
include  mass,  volume,  density,  number  and  in  general  any  scalar 
quantity.  Thus  the  number  of  particles  in  an  arbitrary  volume  or  the 
rate  of  particle  impact  on  a  dectector  moving  at  constant  velocity  are 
independent  of  the  state  of  motion  of  the  observer.  Of  course,  if 
the  state  of  motion  of  the  detector  is  changed,  it  sees  a  different 
particle  distribution  in  velocity  space  as  discussed  above.  In  general 
it  is  easier  to  compute  impacts  on  a  detector  at  rest  than  it  is  to 
compute  impacts  when  both  detector  and  particles  are  moving.  Eq  2.15 
permits  the  computation  of  the  particle  distribution  in  the  coordinate 
system  for  which  the  detector  is  at  rest.  In  this  coordinate  system, 
the  definitions  of  Table  1  are  made. 


2.17 

2.18 
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Table  1.  Definition  of  Various  Physical  Quantities  in  Terms 
of  the  Density  Function  in  Six-Space 


1.  Total  number  of  particles 


N (r, v)drdv 


N 


2.  Density  of  particles 
(Number/m^) 


N(r,v)  dv 


=  N(r) 


3.  Counting  rate  by  a  unit 
isotropic  detector 


4.  Counting  rate  by  a  unit 
Surface  detector  with 
normal  n. 


Ill 

III 


N(r, v)vdv 


N(r,v)v*n 


=  <t>(r) 

dv  =  <Kr,n) 


Because  of  its  association  with  idealized  particle  distributions 
with  respect  to  r  and  v,  it  is  useful  to  enumerate  some  of  the 
properties  of  the  Dirac  6  function  and  to  use  it  to  represent  several 
distributions  of  common  interest.  A  graphical  depiction  of  the  6 
function  and  others  related  to  it  are  shown  in  Figure  2.  By  definition, 
if  the  range  of  integration  includes  the  origin. 


x=a 


J  f(x)  6(x)dx  =  f(0).  2.20 

x=a 

It  follows  that 


r  *=b 

1  6(x)dx  =  1 

2.21 

If 

^  x=a 

x  =  y-a 

2.22 

so  that 

y  =  x  +  a 

2.23 

then 

y=c  x=c-a 

J  f(y)  >S(y-a)dy  =  J  f(x+a)  6(x)dx  =  f(a) 

2.24 

y=b  x=b-a 
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Figure  2.  Some  special  functions.  The  derivative  of  function  1  is 
function  2  etc.  Function  3  is  the  Dirac  6  function. 
Functions  3  and  4  are  shown  as  being  finite  for  graphical 
purposes,  but  actually  extend  infinitely  far. 


Thus,  since  the  variable  of  integration  is  immaterial. 


,  x=c 

L 

f(x)  6(x-a)dx  =  f(a). 

2.25 

Suppose  that 

y  =  F (x) 

2.26 

and 

x  =  G(y). 

2.27 

Then 

x=c 

r 

y-C 

J  6[x-a]dx  = 

x=b 

J  6  [G(y)  -  G(A)](dG/dy)dy  =  1. 

y=B 

2.28 

where  y=A  when  x=a. 

Thus 

<5(y-A)dy  = 

6[G(y)  -  a](dG/dy)dy 

2.29 

or 

6(y-A)  = 

6[G(y)  -  a]dG/dy. 

2.30 

As  an  example,  suppose 

that 

2 

V  =  X 

2.31 

and 

1/2 

x  =  y 

2.32 

and  that  we  wish  to  evaluate  the  integral 

P 

6(x-a)dx  =  a3. 

2.33 

With  a  change  in  variables 
integral  becomes 

according  to  Eqs  2.31  and  2.32,  the 

J>2 

<5[y1/2-a][(l/2)y’1/2]dy. 

2.34 

2 

By  using  Eq  2.30  with  the  definition  A  =  a  ,  this  becomes 

/>2 

6(y-A)dy  =  (A) 3/2  =  a3. 

2.35 
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The  6  function  concept  can  be  extended  with  care  to  more  than 
one  dimension.  For  example,  the  definition, 


«(r)  =  6(x)  6(y)  6(z), 


permits  the  expressions 


III 


F(r)  6(r)  dx  =  F(0), 


and 


F (r)  6 (r-rQ)  dx  =  F(rQ)  . 


2.36 


2.37 

2.38 


It  should  be  noted  that  <Sr  defined  in  spherical  coordinates  is  not 
equivalent  to  6(x)  6(y)  6 ( z)  although  there  are  strong  similarities. 

As  a  rule,  it  is  easier  to  deal  with  continuous  rather  than  discrete 
distributions  because  of  the  difficulty  of  mapping  the  6  function  from 
one  coordinate  system  to  another,  especially  if  there  are  several 
variables. 


3.  TECHNIQUES  FOR  DEALING  WITH  MONODI RECTI ONAL  FLOWS 

In  dealing  with  particle  fields  which  are  locally  monodirectional , 
i.e.  all  particles  in  the  neighborhood  of  a  point  in  r  space  are 
traveling  in  the  same  direction,  the  Liouville  theorem  and  the  Dirac 
delta  function  are  not  particularly  useful.  The  photons  from  a 
distant  star,  a  collection  of  meteoroids  moving  in  a  loosely  connected 
cloud  through  the  solar  system,  or  the  droplets  from  a  sprinkler 
system  are  examples  of  locally  monodirectional  flows.  Suppose  as 
shown  in  Figure  3  that  a  point  source  emits  particles  at  a  constant 
rate  S(a,B),  measured  in  particles/ (steradian  second),  into  the  element 
of  solid  angle  sina  dadB.  If  the  flow  is  time  independent,  at  any 
point  on  a  particle  trajectory  we  must  have 


s(a,3)  sinadadB  =  J  dA  3.1 

2 

where  J  is  the  current  density,  measured  in  particles/m  sec,  and  dA 
is  an  element  of  area  normal  to  the  particle  trajectories.  Solving 
for  J  we  find 


J  =  S(a,B)  sinadadB/ dA  .  3.2 

The  central  problem  in  monodirectional  flows  is  to  find  the  magnitude 
and  direction  of  the  current  J  in  a  convenient  coordinate  system. 
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In  Figure  3,  three  coordinate  systems  may  be  visualized.  The 
source  location  at  P  may  be  viewed  as  a  coordinate  origin  from  which  it 
is  convenient  to  describe  particle  emission  in  terms  of  S(a,B).  An 
observer  at  point  0  may  wish  to  measure  J(r)  in  terms  of  convenient 
spherical  or  Cartesian  systems  with  origin_at  point  0.  For  convenience 
in  calculating  the  vector  current  density  J,  it  would  be  convenient 
to  have  a  coordinate  system  based  on  the  normal  and  tangent  unit  vectors 
of  particle  trajectories. 

Instead  of  the  point  source  of  Figure  3,  a  parallel  stream  passing 
an  attractive  center  is  shown  in  Figure  4.  Following  the  same  line  of 
reasoning  as  before, 

J(r)  =  J(p, 6)  a  da  dS/dA  3.3 

where  J(p,6)  is  the  current  density  distribution  at  a  point  distant 
from  the  attractive  center  and  J(r)  is  the  perturbed  current  density 
near  it. 

In  general  it  is  easy  to  establish  an  element  of  area  dA'  which 
contains  the  point  located  by  r,  but  it  must  be  projected  onto  the 
plane  normal  to  the  trajectory  such  that 

dA  =  dA'  cose,  3.4 

where  e  is  the  angle  between  the  normals  of  dA  and  dA'.  In  spherical 
coordinates,  for  example,  any  one  of  three  elemental  areas  can  be 
chosen,  videlicet, 


dA'  =  r  dr  d0  <p  , 

3.5 

or 

dA'  =  r  sin©  dr  dip  8  , 

3.6 

or 

2 

dA  =  r  sinG  do  d<j>  r  , 

3.7 

where  r  0  and  £  are  the  appropriate  unit  vectors  defined  by  the 
spherical  coordinate  system  to  be  normal  to  the  elemental  areas. 


The  trajectory  or  orbit  of  a  particle  may  be  expressed,  in  terms 
of  the  variables  defined  by  Figure  3,  as  the  intersection  of  two 
surfaces, 

F(r,9,<f>,a)  =  0  3.8 

and 

G(r,  0,  <j>,  8)  =  0.  3.9 
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A  depiction  of  monodirectional  flow  from  a  point  source 
located  at  the  point  P.  The  areas  dA  include  the  trajectories 
of  particles  emitted  into  the  solid  angle  sina  dadg. 


Figure  4.  A  depiction  of  a  parallel  stream  interacting  with  an 

attractive  center.  The  variable  "a"  may  be  defined  as  the 
impact  parameter. 
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A  third  coordinate  surface  may  be  defined  by  the  distance  the  particle 
has  traveled  from  the  source.  Since  the  functions  VF  and  VG  are 
normal  to  the  trajectory,  they  can  be  used  to  define  a  unit  vector  t 
which  is  tangent  to  the  trajectory  by  the  equation, 

t  =  (VF  x  VG)/ | VF  x  VG|.  3.10 


If  0  and  <|>  are  chosen  to  define  the  elemental  area  dA'  of  Eq  3.4, 

2  -  - 

dA  =  dA'  cose  =  r  sin0  d 0d d>  r*t,  3.11 

and  Eq  3.2  becomes 

_  2  -  - 

J(r)  =  s(a,B)  sinadadB/[r  sin0  d0dj>  r*t],  3.12 

where 

cos  e  e  r-t.  3.13 

Likewise,  Eq  3.3  could  be  written  as 

J(r)  =  J(p ,  6)  a  da  d6/ [r2  sin0  d0  d<t>  r*t]  .  3.14 


In  general,  if  equations  can  be  found  connecting  a  and  B  of  Eq  3.12 
to  6  and  <J>,  we  can  write 


dadB 


3a/ 30 
3a/  3d> 


36/30 
36/  3d> 


d0d<}> 


3. 15 


E  D(||)  d 6d4>  3.16 

where  the  determinant  D  is  the  Jacobian  transformation  connecting 
coordinate  differentials.  This  permits  Eq  3.12  to  be  written  finally 

as 

J(r)  =  s(a,6)  sina  D(^)/[r2  sin0  ?•£].  3.17 

Eq  3.17  contains  formally  the  magnitude  J(r)  and  direction  t  of  the  vector 
current  density  at  an  arbitrary  point  in  space  with  t  given  by  Eq  3.10. 
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If  particles  are  neither  created  or  destroyed  during  steady-state 
flow,  the  continuity  equation  may  be  written  as 

V-J  +  =  0  3.18 

and  since  by  definition  of  steady-state  flow 

3N/3t  =0  3.19 

we  must  have 

V-J  =  0  3.20 

or  divergenceless  flow.  The  preceding  discussion  of  course  is  based 
on  the  assumption  of  divergenceless  flow  and  a  corollary  that  the  flow 
contained  in  a  bundle  of  streamlines  is  the  same  across  any  area 
defined  by  the  streamlines  and  normal  to  them.  Thus 


JO^Afrp  =  J(r2)A(r2)  3.21 

or 

NCFpvcTpAffp  =  N(F2)V(>2)A(>2)  3.22 

for  any  two  points  on  a  particle  trajectory. 

4.  UNIFORM  FORCE  FIELDS 

To  illustrate  the  usefulness  of  the  techniques  developed  in  the 
foregoing  pages,  two  standard  problems,  dealing  with  idealized  waterfalls 
and  sprinklers,  will  be  examined  from  several  viewpoints. 

We  will  assume  in  the  waterfall  problem  that  we  have  a  vertical 
fall,  in  a  uniform  gravilational  field,  of  particles  which  do  not 
interfere  with  each  other  and  that  the  flow  is  constant  in  time. 

If-Z  is  the  distance  measured  downward  from  the  coordinate  origin, 

Eq  3.22  yields  the  solution, 

N(Z)v(Z)  =  N(0)v(0)  4.1 

where  the  conservation  of  energy  equation, 

v2(Z)/2  +  gZ  =  v2(0)/2,  4.2 

yields  v(Z)  as 

v(Z)  =  [-2gZ  +  v2(0)]1/2.  4.3 
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In  these  equations  v(0)  is  the  speed  when  Z  is  zero  and  g  is  the 
gravitational  acceleration. 

The  waterfall  problem  could  be  approached  more  generally  by 
assuming  a  distribution  of  velocities  given  by 

N(0,v)  =  K  ,  v1<v<v2.  4.4 

The  Liouville  theorem  permits  N(Z,v)  to  be  written  as 

N(Z,v)  =  K  ;  v1(Z)<v<v2fZ)>  4.5 

where  from  Eq  4.3, 

Vj(Z)  =  [vj  -  2gZ] 1/2  4.6 

and 

v2CZ)  =  [v2  -  2gZ] 1/2  4.7 

By  definition, 

r  V2^Z^ 

N(Z)  =  I  N(Z,v)dv 

Vj  (Z) 

=  K  [v2C2)  -  Vj(Z)].  4.8 

In  the  limit  as  v^  approaches  v 2,  v^  and  v2  can  be  replaced  by  v(0) 
and  v2  -  Vj  by  A  v(0) .  From  Eq  4.8 

N(Z)  =  KAv(Z)  4.9 

and  from  Eqs  4 . 6  and  4 . 7 

Av(Z)  =  v(0)  Av(0)/[v2(0)  -  2gZ]1/2  4.10 

so  that 

N(Z)  =  K  v(0)  Av(0)/[v(0)2  -  2gZ]1/2.  4.11 

From  Eqs  4 . 1  and  4 . 3 

N (Z)  =  Nov(o)/[v(0)2  -  2gZ] 1/2  4.12 

If  K  v(0)  Av(0)  is  identified  with  NfO)v(O),  Eqs  4.11  and  4.12  ate  the 
same,  thus  demonstrating  the  value  of  the  Liouville  theorem. 
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In  the  more  complex  sprinkler  problem,  let  us  assume  that  a  water 
sprinkler  emits  droplets  according  to  the  density  function  S(v,a,(3). 
This  means  that  the  expression  S(v,a,3)  v^  dv  sinadadg  defines  in 
velocity  space  the  number  of  particles  emitted  each  second  having 
speeds  between  v  and  v  +  dv  and  direction  lying  in  the  solid  angle 
sinadadg.  The  problem  is  to  find  the  rate  at  which  particles  fall  on 
an  arbitrary  unit  area  in  the  plane  containing  the  sprinkler.  The 
coordinates  of  the  problem  are  defined  in  Figure  5.  The  choice  of 
the  coordinate  system  permits  the  equality  of  6  and  <f>. 

Starting  with  Newton’s  second  law,  the  trajectory  for  a  particle 
emitted  with  speed  v  in  the  direction  defined  by  a  and  6  is  given  by 

< p  =  $  =  constant  4.12 

P  =  v  sinu  t  4.14 

2 

Z  =  v  cosa  t  -  gt  /2  4.  IS 

where  t  is  the  time  elapsed  after  launch  and  g  is  the  gravitational 
acceleration. 

The  particle  hits  the  surface  when  Z  equals  zero  so  that 

2  2 

p  =  (2v  /g)  sina  cosa  =  (v  /g)  sin  2a  4.16 


If  it  is  assumed  that  the  emission  at  speed  v  in  the  solid  angle 
sinadad3  falls  into  the  elemental  area  pdpdiji  and  it  is  recalled  that  $ 
and  3  are  identical, 

S(a,<J>)  sinadad<f>  =  J  cos  e  pdpd<j>  4.1" 


From  Eq  4.16, 


so  that 


T  S(a,4>)  sinada  _  T 

J  cos  e  =  — — — - t—  =  J7 

p  dp  Z 


=  2(v2/g)  cos  2a 


j  _  S(a,()>)  sina _ 

2  2(vVg2)  sin  2a  cos  2a 


S(a,  <k)  sina 


(v4/g2)  sin  4a 


Li 


where  from  Eq  4.16 


a  =  (1/2)  sin  *  (pg/v2) ,  0<a£  tt/2  4.21 

For  the  particular  case  of  the  isotropic  source, 

S  =  Q/2tt  4.22 

where  Q  is  the  total  water  flow, 

=  (Q/2tt)  sin  a/[(v4/g2)  sin  4a].  4.23 

For  this  particluar  case  there  is  symmetry  in  the  trajectory  and  t  is 
obtained  by  inspection  to  be  given  by  equation, 

t  =  -£  cosa  +  i  sina  cosif>  +  j  sina  simj>,  4.24 

which  is  just  a  reflection  of  the  original  direction  of  emission  in  the 

x  y  plane. 

It  should  be  noted  that  dp/da  changes  sign  when  a  =  ir/4,  the  angle 
which  gives  maximum  range.  Thus  each  wetted  area  receives  two  flows, 
one  for  0<a£  tt/4  and  one  for  ir/4  <ol<_  tt/ 2 .  These  are  plotted  in  Figure 
6.  If  it  is  desired  to  produce  a  sprinkler  which  wets  uniformly, 

Eq  4.23  tells  us  that  we  must  have 

S (a,  $)  =  k[(v4/g2)  sin  4a]/sina,  4.25 

where  k  is  a  constant  of  proportionality. 

5.  CENTRAL  FORCE  FIELDS  AND  ISOTROPIC  DISTRIBUTIONS 

When  manned  space  flight  was  in  its  early  stages  there  was  a 
considerable  interest  in  the  distribution  of  meteoroids  about  the  earth. 

If  the  meteoroid  distribution  about  the  earth  were  assumed  to  be 
isotropic  in  direction,  a  theoretical  treatment  by  S.  F.  Singer1  yielded 
the  density  of  meteorids  as  a  function  of  velocity  at  infinity  and 
distance  from  the  earth.  If  the  meteoroids  were  thought  to  consist  of 
a  locally  monodirectional  stream  flowing  by  the  earth,  a  theoretical 
treatment  by  R.  D.  Shelton^  et  al  gave  the  density  and  flux  distributions. 


^S.  F.  Singer,  Nature,  Oat  78,  1961,  pp  321-323. 

2 

?.  D.  Shelton,  et  al.  Astronomical  Journal  70  No. 
170. 


2,  Mar  1965,  pp  166- 
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which  sprays  isotropically.  Each  wetted  point  receives 
from  two  directions. 


The  following  discussion  of  isotropic  meteoroid  distributions 
differs  slightly  from  that  of  Singer,  but  illustrates  the  use  of 
techniques  previously  developed  and  obtains  the  same  results.  Let 
the  meteoroid  population  far  from  the  earth  be  given  by 

S(r)  =  K  =  S(”)  ;  Vj-cvcv^  5.1 

The  density  at  infinity  is  given  by  definition  as. 


N(r)  =JJJ 


K  v2  dv  sine  de  d<J>  =  C4 ttK/3)  CV^  -  v^)  =  N (<*>).  5.2 


The  flux  measured  by  an  isotropic  unit  detector  would  be  given  as 


r  r  r 

=  K 
JJ  J 


3  4  4 

v  dv  sine  ded<(i  =  ttK(v2  -  v^)  =  <t>(°°) 


In  the  above  formulae,  v,  6  and  <j>  are  coordinates  in  velocity  space. 

The  conservation  of  energy  for  particles  passing  near  the  earth  requires 
that 

v2  =  v2(r)  -  2GM/r,  5.4 


where  v^  is  the  particle  speed  at  infinity,  v(r)  is  the  particle  speed 

after  it  has  moved  to  a  location  r,  G  is  the  gravitational  constant 
and  M  is  the  mass  of  the  earth.  The  Liouville  theorem  and  Eq  5.1 
require  that 

S(r)  =  K  ;  Vl(r)<  v<  v2(r)  5.5 

where  from  Eqs  5.1  and  5.4 

vx(r)  =  [v2  +  2GM/r]1/2  5.6 

and 

v2(r)  =  [v2  +  2GM/r]1/2.  5.7 
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Therefore 


Nfr) 


dv  sine  d0d(f> 


(4ttK/3)  {  [v^  +  2GM/r]3/2  - 


5.8 


[Vj  +  2GM/r] 3/2 } 
and 

<Kr)  =  ttK  {[v2  +  2GM/r]2  -  [v2  +  2GM/r]2}  5.9 

In  the  limit  as  Vj  and  v2  approach  (i.e.,  the  distribution  is 

monoenergetic  and  isotropic),  and  their  difference  becomes  AV,  Eqs  5.8 
and  5.9  may  be  written  as 


N(r)/N(°°)  =  vfrj/v^  5.10 

and 

<t>(r)/<K°°)  =  [v(r)/vj2  5.11 

where 

v(r)  *  [v2  +  2GM/r]1/2.  5.12 

The  meteoroid  distribution  is  diminished  near  the  earth  because 
the  earth  stops  the  motion  of  meteoroids  whose  trajectories 
intersect  it.  By  using  the  principle  of  conservation  of  angular 
momentum,  as  shown  in  Figure  7,  we  can  write 

R  v(R)  =  r  v(r)  sin9m  S.13 

so  that 

Sin0m  =  R  v(R)/r  v(r).  5.14 

The  fractional  solid  angle  f  available  at  point  P  is  given  by 


1  r$=2  TT  ^0=7 

47  J  J 

J4>=0  0=0 


■  0=  TT-  0 


m 


sine  d9d$. 


5.15 
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1 


Figure  7.  The  geometry  of  grazing  orbits.  The  conservation  of  angular 

momentum  for  a  trajectory  which  grazes  the  earth  requires 

that  r  v(r)  sinO  =  R  V(R) . 
m 


r 


From  Eq  5.13 

cos  0m  =  [l-sin20m]1/2  =  [l-R2v2(R)/r2v2(r)]1/2  5. *6 

so  that,  using  Eqs  5.4,  5.15  and  5.16, 

f  =  1/2  {1  +  [1- (R/r) 2 (v2  +  2GM/R) / (v2  +  2GM/r)]1/2}.  5.17 

As  expected,  f  goes  to  1/2  near  the  earth  and  to  unity  at  large 
distances. 

Using  Eqs  5.10,  5.11  and  5.17  we  can  write 


N(r)/N(“)  =  f  vfrj/v^  5.18 

4>(r)/<K»)  =  f[v(r)/vj2  5.19 

It  is  convenient  to  modify  the  preceding  equations  by  the  definitions 

R  =  1  5.20 

vg  =  2GM/R  5.21 

and 

v  =  vjve  5.22 


so  that  r  is  measured  in  earth  radii  and  v  is  measured  in  terms  of 

00 

the  escape  velocity.  Eq  5.4  now  becomes 

v2(r)  =  v2  +  (2GM/R) (R/r)  =  v2  +  1/r  5.23 

and  Eqs  5.10  and  5.11  becomes 

N(7)/N(»)  =  [1  +  1/ (rV2) ] 1/2  (1/2)  {1  +  [l-(l/r)2(V2  +  1)/(V2  +  1/r) ] 1/2 }5. 
$(7)/<K»)  =  [1  +  1/ (rV2)  ]  (1/2)  (1  +[l-(l/r)2(V2  ♦  1)/(V2  +  1/r) ] 1/2}5 . 25 

Note  that,  as  r  goes  to  infinity,  these  ratios  go  to  unity  as  expected 
physically.  As  V  becomes  large  and  the  earth  has  no  influence,  these 
ratios  are  diminished  by  the  fraction  of  the  total  solid  angle  subtended 

2  2 

by  the  earth  at  distance  r,  namely  R  /r  .  Eq  5.23  is  the  same  as 
Singer's  Eq  8,  except  for  a  difference  in  the  definition  of  variables. 

Eqs  5.24  and  5.25  have  been  computed  and  plotted  in  Figures  8  and  9. 


There  have  been  speculations  about  populating  various  orbits  with 
fragments  as  an  anti-satellite  measure.  As  a  crude  example,  suppose 
that  it  is  desired  to  put  fragments  into  orbits  from  100-300  miles 
above  the  earth  in  the  zone  extending  to  20°  latitudes  north  and  south 
of  the  equator.  The  volume  of  this  region  is  given  by  the  formula 

3  3  20  3 

V  =  (4ir/3)  (r^  -  rp  cos  0Q  =  1.7  x  10  M  ,  5.26 


where  r0  is  the  outer  radius,  r,  is  the  inner  radius  and  0  is  the 
2  1  °  6 
angular  extent  of  the  zone  measured  from  the  equator.  If  10  particles 
are  injected  into  a  retrograde  orbit  so  that  they  encounter  the  targets 
with  maximum  speed  and  flux,  the  speed  of  impact  is  roughly  twice  the 
orbital  velocity  or  15  km/sec.  The  impact  rate  is  computed  to  be 

<{>  =  N  v  =  CIO6/ 1 . 7  x  1020)  x  16  x  103  =  9.4  x  10~U  p/m2  sec  5.27 

or  one  impact  on  a  square  meter  of  surface  in  every  337  years.  Thus  it 

is  seen  that  a  cloud  of  projectiles  is  not  an  effective  way  to  kill 

hostile  satellites  and  that  a  chance  encounter  with  debris  left  in 
space  is  vanishingly  small,  being  about  one  in  1,000,000  years  if  we 
assume  that  300  pieces  are  still  in  orbit. 

A  necessary  condition  for  a  particle  to  be  injected  into  an  orbit  which 
misses  the  earth  is,  from  Eq  5.13,  that 

R  v(R)  £  Rt  v(Rp  sin  @m  5.28 

From  the  requirement  that  energy  be  conserved  (Eqs  5.4  and  5.21), 

v2(R)  -  v2  =  v2  (Rp  -  v2  (R/Rp. 

From  Eqs  5.28  and  5.29, 

v2(Rp  >  [(R2/Rp(RrR)/R2  sin2  0m  -  R2)]v2, 

which,  when  0  equals  ir/2,  becomes 
m 

v2  (Rp  1  [(R2/Ri)/(Ri  +  R)v2  .  5.31 


5.29 


5.30 
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5.32 


Thus  we  may  write 

v2  (R2/R.)/(R.  +  R)  1v2(Rj)  <  v2  (R/R.)  . 

The  quantity  on  the  right  is  the  escape  velocity  for  a  particle  located 
a  distance  R.  from  the  earth.  At  5  earth  radii,  for  example,  Eq  5.31 
becomes 


ve//30  <  v  (R.)  1  ve//5  . 


5.33 


6.  MONOENERGETIC  MONODIRECTIONAL  FLOWS  PAST  A  GRAVITATIONAL  CENTER 

If  the  earth  should  encounter  a  large  loosely  connected  cloud  of 
dust  and  meteoroids,  it  would  appear  to  an  observer  on  the  earth  that 
a  stream  of  particles,  with  parallel  trajectories  at  large  distances 
on  the  approach  to  the  earth,  was  flowing  past  the  earth  and  being 
influenced  by  its  gravitational  field,  as  shown  in  Figure  10.  Each 
particle  orbit  is  characterized  by  a  velocity  v^  at  large  distances  from 

the  earth,  an  impact  parameter  "a",  and  an  azimuthal  angle  <t>.  The 
impact  parameter  associated  with  a  particle  may  be  defined  as  the 
distance  of  closest  approach  of  the  particle  to  the  center  of  the  earth 
if  the  gravitational  field  were  absent  and  the  particle  traveled  in  a 
straight  line.  A  particle  trajectory  can  be  defined  by  the  equations, ^ ^ 

F  i  r  [1  +  e  cos  (0-9')]  -  B  =  0  6.1 

and  G  =  ip -<p'  =  0,  6.2 

where 

B  i  2  a2  v2/R  v2  6.3 

oo  e 


and 

E  =  [1  ♦  (4  a  vJ/(R  ve)]  .  6.4 

and  vg  is  the  escape  velocity  and  R  is  the  radius  of  the  earth. 

The  constant  6'  is  defined  by  the  requirement  that  r  be  infinite  when 

0  takes  on  the  values  of  zero  and  0  .  From  Eq  6.1  it  follows  that 

m  M 

cos  0'  =  (-1/e)  =  cos  (0  -0')  6.5 


George  Joos,  Theoretical  Physics,  Hafner  Publishing  Co.,  NY,  1950,  p  89. 

Herbert  Goldstein,  Classical  Mechanics,  Addiaion  Wesley  Publishing  Co. , 
Inc.,  Reading,  MA,  1959,  p  77. 
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and  that 


0  / 2  =  e'  =  .-os'  (-1/e). 


If  the  definitions  are  made  that 


and 


R  =  1 


U  =  v  /  v  , 

<*r  e 


Equations  6.3  and  6.4  can  be  written  as 

2  2 

B  =  2  a  U 


and 


M  .  2..4.1/2 

e  =  [1  +  4  a  U  ] 


6.6 


6.7 


6.8 


6.9 


6.10 


The  calculation  of  the  vector  current  density  J(r,9)  can  proceed 
directly  from  the  formalism  leading  to  Eqs  3.10  and  3.14.  If,  for 
example,  we  identify  the  same  particle  flux  through  an  element  of  area 
a  da  d<j>  at  large  distances  to  the  left  of  Figure  10  and  through  an  element 
2 

of  area  r  sin9  d0  d4>  located  at  an  arbitrary  point  fr,  0,<j>)  we  can  write 
the  equation, 

a  da  d$  =  J(r,9)  r2  sine  de  d$  (r*t).  6.11 

Since  the  Jacobian  connecting  the  differential  areas  reduces  to  3a/36, 
we  can  write 


J(r,0)  =  [(J^  a)/(r2  sin6  •  r*t) ] (3a/36) .  6.12 

^  A  A  2 

The  vector  scalar  product  r*t  serves  to  project  the  area  r(r  sin9  de  d<j>) 
onto  a  plane  normal  to  the  trajectory,  the  direction  of  which  is  defined 
by  t.  In  a  like  manner,  by  choosing  a  different  element  of  area 
9  r  sine  dr  d<{>,  we  can  write, 

J(r,0)  =  [(Jw  a)/(r  sin9  •  9*t) ] (3a/3r) .  6.13 

By  combining  Eqs  6.1,  6.S,  6.9  and  6.10,  the  values  of  "a"  associated 
with  a  given  r  and  0  can  be  found.  Using  the  trigonometric  identity, 

cos(9-9')  =  cos6  cos0'  +  sin9  sin9'  6.14 

we  can  use  Eq  6.5  to  write 

2  1/2 

e  cos(e-e')  =  -COS9  +  (e  -1)  sine,  6.15 


39 


Figure  10.  Flow  of  a  parallel  stream  of  meteoroids  from  left  to 
right  past  the  earth.  Points  (1)  shielded  by  the 
earth  will  see  no  meteoroids.  Some  points  (2)  will 
receive  from  one  direction,  and  other  (3)  will  receive 
from  two  directions.  Points  along  the  symmetry  axis  (4) 
will  receive  from  all  directions  lying  on  the  cone 
defined  by  the  angle  5. 
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where  we  have  used  the  trigonometric  identity, 

sinO'  =  (l-cosV)1/2  =  (e2-l)1/2/e.  6.16 

Thus  Eq  6.1,  with  the  help  of  Eqs  6.9,  6.10  and  6.15  can  be  written  as 

2  1/2  22 

r{l  +  [-cos0  +  (e  -  1)  '  sin0]  -  2  a  U  }  =  0  6.17 

or 

r(l  +  [-cose  +  2a  U2  sine]  -  2  a2  U2 }  =0.  6.18 

This  equation  can  be  solved  for  "a"  to  yield  the  result, 

a  =  C 1  / 2)  {  r  sin6  +_  [r2  sin2  6  +  2r(l-cos6)/(J2]1^2}.  6.19 

If  follows  that 

3a/3r  =  (1/2)  {sin6  +  [rsin^  6  +  (l-cos6)/U2]  }  6.20 

[r2sin2  6  +  2r(l-cos6)/U2] 1//2 

and 

„  2  2 

8a/3e  =  (1/2)  {rcose  +  [r  sine  cos6  +  (r  sin6)/U  ]  ).  6.21 

22  21/2 

[r  sin  6  +  2r(l-cos6)/U  ] 

By  using  our  previous  definition  of  t  as 

t  =  VF  x  VG/|VF  x  VG|  6.22 

Eqs  1  and  2  yield  the  result  in  spherical  coordinates  that 

VF  =  {-r/r2  +  (6/Br)e  sin(6-6')  6.23 

and 

VG  =  {[l/(r  sin6)]4>,  6.24 

so  that 

t  =  [6/r  +  (r/B) e  sin(e-6 ') ]/ [1/r  +  (e/B)  sin  (6-6')]  '  .  6.25 

Thus 

r-t  =  [(e/B)sin(6-6')]/[l/r2  +  (e/B) 2sin2 (6-6 ') ] 1/2  6.26 
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1 


and 


0-t  =  [(l/r)]/[l/r2  +  (e/B)2sin2(6-6')]1/2.  6.27 

where  from  Eqs  6.9  and  6.10 

e/B  =  (1  +  4  a2U4)1/2/(2  aV) .  6.28 

Figure  11  illustrates  the  relationship  of  two  trajectories  which 

intersect  at  r  and  0.  The  particle  flux  <j>(r,0)  at  r  and  0  results  from 

the  sum  of  two  flows,  one  associated  with  a  and  the  other  with  a-,  so 

’  +  ’ 

that  we  can  write 


J (r,  0)  =  |  J+(r,  0)  |  +  |  J_  (r, 2ir-0)  |  6.29 

where  J+(r,0),  for  example,  is  calculated  from  Eq  6.12  using  the  values 

of  "a"  and  3a/ 30  associated  with  the  positive  signs  of  Eqs  6.19  and  6.21 
and  the  same  value  of  "a"  to  compute  e/B  of  Eq  28  and  r*t  in  Eq  6.26. 

Because  of  shielding  by  the  earth,  J-(r, 0)  or  both  J+(r,9)  and 

J-(r,0)  may  not  reach  a  point.  Trajectories  with  impact  parameters  less 
than  that  associated  with  a  grazing  trajectory  will  intersect  the  earth 
and  fail  to  contribute  to  the  particle  flux  thereafter.  If 


0  <  e' 


and 


r  >  1 

the  trajectory  has  not  reached  its  closest  approach  to  the  earth. 

0  >  0' 


6.30 

6.31 
If 

6.32 


and 

rc  <  1  6.33 

where  rc  is  the  point  of  closest  approach,  the  flux  contribution  from 
that  trajectory  is  zero. 


The  distance  of  closest  approach  rc  is  easily  calculated  from 
the  conservation  of  angular  momentum  and  energy, 

U  a  =  U(l)  6.34 

g 
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and 


U2  =  U2(l)-1, 


6.35 


where  the  definition  of  Eqs  6.7  and  6.8  are  used,  U(l)  is  the  speed 
at  the  grazing  point,  and  a  is  the  impact  parameter  associated  with 

a  grazing  trajectory.  Eliminating  U(l)  and  solving  for  ag  we  have  the 

result. 


ag  =  (U2  +  1) 1/2/U. 


6.36 


Thus  any  trajectory  with  "a"  less  than 


"a  "  will  intersect  the  earth, 
g 


The  tests  and  procedures  for  computing  <|>(r,0)  are  contained  in 
Appendix  1.  The  isoflux  plots  for  <t>(r,0)  are  shown  in  Figures  12-14 
for  values  of  U  equal  to  0.5,  1  and  2.  The  infinite  flux  along  the 


symmetry  axis  implies  that,  during  a  heavy  meteor  shower,  it  might  be 


possible,  under  the  right  lighting  conditions,  to  see  the  symmetry 
axis  as  a  fluctuating  light  spot  in  the  sky  on  the  opposite  side  of  the 


earth  from  the  direction  of  meteoroid  arrival. 


7.  FOCUSING  OF  PHOTONS  BY  A  SCHWARZSCHILD  FIELD 


The  preceding  methods  lend  themselves  readily  to  the  treatment  of 
stellar  phenomena  associated  with  the  bending  of  light  rays  by 
gravitational  fields.  The  light  rays  are  viewed  as  a  swarm  of  photons, 
each  following  its  individual  trajectory  as  it  moves  through  the  spherically 
symmetric  field,  first  computed  by  Schwartzschild,  created  by  a  massive 
body.  According  to  the  general  theory  of  relativity^,  light  rays  passing 
by  a  massive  celestial  body  are  bent  by  an  amount  ip  according  to  the 
formula, 

ip  =  4GM/ac2H  K/a  7.1 

where  G  is  the  universal  gravitational  constant,  M  is  the  mass  of  the 
body,  usually  considered  to  be  a  star,  c  is  the  velocity  of  light,  and 
"a"  is  the  impact  parameter.  It  is  conjectured,  although  it  has  never 
been  observed,  that  an  observer,  in  line  with  two  stars  as  shown  in 
Figure  15,  will  see  the  distant  star  magnified  greatly  by  the  lense 
effect  of  the  Schwartzschild  field  associated  with  the  nearer  star. 


5P.  G.  Bergmann,  Introduction  to  the  Theory  of  Relativity,  Prentice  Hall 
Inc.,  Englewood  Cliffs,  NJ,  1942,  p  221. 
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Figure  12.  An  isoflux  plot  for  a  parallel  stream  of  meteoroids  moving 
by  the  earth  at  half  the  escope  velocity  (U  =  1/2).  The 
dotted  lines  are  grazing  orbits,  which  demarcate  the  region 
A  receiving  no  flux,  the  region  B  receiving  flux  in  one 
direction,  and  the  region  C  receiving  flux  in  two  directions. 
Intensities  are  calculated  by  assuming  unit  flux  at  infinity . 
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Figure  13.  An  isoflux  plot  for  U  =  1.  For  an  explanation  of  the  symbols, 
see  Figure  12. 


Light  leaving  the  source  at  (1)  in  the  direction  in  Figure  15 
will  travel  in  a  direction  defined  by  a  according  to  the  equation/ 


r  cosfe-ej  ')  =  a  7.2 

where  the  variables  and  parameters  are  as  defined  in  Figure  15.  After 
interaction  with  the  star  located  at  (3),  the  light  will  travel  in  a 
straight  line  according  to  the  equation. 


r  cos (0-02^)  =  a  .  7.3 

From  Eq  7.1  and  Figure  15,  the  amount  of  bending  in  the  light  ray  is 
given  by  the  equation, 

V  =  e2'  -  0i'  =  K/a.  7.4 

Also,  from  Figure  15, 


0j'  =  0/2)  -a,  7.5 

so  that  from  Eqs  7.4  and  7.5, 

@2^  =  'll*  dj  '  =  ip  +  ir/2-a.  7.6 

Using  the  identity  that 


cos(x-tt/2)  =  sin  x  7.7 

and  Eqs  7.5  and  7.6  in  Eqs  7.2  and  7.3,  we  have 


r  sin(0  +  a)  =  a 


and 


r  sin(0-i//  +  a)  =  a. 


7.8 


7.9 


We  will  now  make  the  approximation,  to  be  demonstrated  later, 
that  o,6,i|/,0,  0j  and  n-d2  are  very  small  so  that  Eqs7.8  and  7.9  can 

be  written  as 


rjOj  +  a)  =  a 


7.10 


Guy  Mayer,  Collective  Features  of  Gravitational  Focusing ,  Optics 
Communications ,  Vol.  16,  No.  3,  pp  317-319,  Mar  1976.  This  reference 
takes  an  optical  approach  to  the  same  problem  and  has  an  excellent 
bibliography  on  light  focusing  by  gravitational  fields. 
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Figure  15.  The  lense  effect  associated  with  gravitational  fields.  The 

observer  at  (2)  sees  light  from  a  distant  star  at  (1)  focused 
by  a  nearer  star  at  (3)  so  that  thestar  at  (1)  appears 
magnified  and  distorted.  The  angles  are  greatly  exaggerated. 
The  rays  arriving  at  (2)  from  sources  at  (1)  and  (3)  have  an 
apparent  angular  separation  6. 
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and 


where  the  identity, 


r2(ir-02  +  ’('-a)  =  a 


sin  x  =  sin  (tt-x)  , 


has  been  used  to  obtain  Eq  7.11  from  Eq  7.9.  The  subscripts  1  and  2 
have  been  introduced  to  designate  the  incoming  and  outgoing  rays.  Eqs 
7.4,  7.10  and  7.11  can  be  used  to  solve  for  a  and  a.  From  them, 

a  =  r^  {[(ir-02)  +  0^  _+  ^[(ir-02)  +  Bj]2  +  4K(rx  +  r2)/ Cr1r2)^2 } 

/[  2(rj  +  r2)]  7.1! 


a  =  —  -  0, 

ri  1 


=  -01  +  r2  {[(ir-02)  +  0X] 
/[2(rx  +  r2)] 


02)  +  Qj]  +  4K(rx  +  r2)/(r 


1  r2  1 


3a_  _  +  r2  +  K^-92)  *  _ 7,15 

902  2frl  +  t2)  ~  2(rx  +  r2)  {  [Ctt-02)  +  Sj] 2  +  4K(rj  +  *2)/rir2)  }1/2 


For  our  purposes  we  will  assume  that  r^  and  r2  are  equal  to  r  so  that 
Eqs  7.14  and  7.15  reduce  to  the  equations, 

a  =  -0j  +  (1/4)  (  [(tt-02)  +  ©j]  1  ^[(if-e2)  +  Qj]2  +  SK/r)^1^2}  7.16 


_9a_  =  1  +  1  +  9ll _ 

802  4  1  4  C[C*-02)  +  0X]2  +  8K/r) 1/2 
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From  Figure  15, 


and  from  Eqs  7.3  and  7.6, 


r 2  sin  6  =  a. 


7. 18 


r^  cos  [  -  'l'  —  C 7T/  2 )  +  a]  = 

From  the  identity  of  Eq  7.12,  Eq 


r2  sin[02-^  +  a]  =  a. 

7.19  can  be  written  in  the  form 


7. 19 


r2  sin[ir-e2  +  i|>-oi]  =  a  .  7.20 

Thus,  comparing  Eqs  7.18  and  7.20, 

6  =  ir-02  +  ip-ct .  7.21 


Using  Eqs  7.4  and  7.14, 

6  =  (tt-62)  +  ((/[^(Bj  +  a)  ]  -a,  7.22 

where  a  is  given  by  Eq  7.14  or  7.16. 

The  angle  6  given  by  Eq  7.22  is  the  apparent  angular  separation 
of  light  arriving  from  point  sources  of  light  (stars)  located  at  point 
(1),  (r^,  6j),  and  point  (3),  the  coordinate  origin,  of  Figure  15.  To 

appreciate  the  optics  of  the  gravitational  field,  it  is  useful  to 
calculate  the  distortion  of  the  circular  stellar  disc  and  to  compute 
the  relative  intensity  of  various  parts  of  the  image  from  the  point 
of  view  of  the  observer  at  point  (2)  of  Figure  15. 

To  perform  the  first  task,  the  stellar  image  is  depicted  as  shown 
in  Figure  16.  The  angles  0^  and  <(> ^  associated  with  an  emission  from 

the  point  P  on  the  circumference  of  the  stellar  disc,  are  related  to 
the  angular  displacement  of  the  star  0  and  the  rotation  angle  cj>  by  the 
equations 

cot  ifrj  =  (0q/0)  esc  <t>  -  cot  <(>  7.23 

and 


0j  =  0  sintfi/sinijij 


7.24 
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where  0  is  the  angle  subtended  by  the  radius  of  the  stellar  disc. 

Because  the  trajectories  lie  in  a  plane,  <}>j  and  $2  are  equal,  and  we  can 

use  the  value  of  6j  in  Eqs  7.16  and  7.21  to  compute  two  values  of  6  for 

each  value  of  0^. 

Suppose  in  Figure  15  that  6j  is  zero  and  that  the  star  at  point 

(1)  emits  an  amount  light  S  sina  dad(}>  into  the  element  of  solid  angle 
sina  dad<(>.  Assume  that  this  light  arrives  at  point  (2)  so  that  an 

—  „  »  2 

amount  of  light  J(r2)  n,t  r2  s^n  92^02^^  Passes  throught  the  element 
2 

of  area  r2  sin  @2d 02d4> •  Equating  these  two  expressions  we  have 

J(r2)/S  =  (sin  a/n-t  r2  sine2) (3a/302)  7.25 

We  can  assume  that  a  and  ir-@2  are  very  small  and  that  n*t  is  unity  so 
that  we  have 

J(r2)/S  =  |{a/[r2(Tr-62)]}  (3a/302) |  .  7.26 

The  radiation  directly  from  the  star,  using  the  inverse  square  law, 
would  be  given  by  the  equation, 

J(r2)/S  =  l/[(2r)2]  =  l/(4r2)  7.27 

assuming  that  r^  and  r2  are  equal  to  r. 

We  could  rotate  the  diagram  of  Figure  15  by  an  angle  -0^  about 

£oinjt  (3)  without  changing  the  physical  relationships  involving 
Tj,  r2,  a,  and  i|>.  Thus  Eq  7.26  should  be  applicable  to  the  more 

general  condition  of  Figure  15  if  we  replace  tt-62  everywhere  by 

[Cir-0  )  +  0^].  As  expected,  "a"  calculated  by  Eq  7.13  is  unchanged 

by  this  rotation.  Although  a  is  changed  by  the  rotation  (-0^)  according 

to  Eq  7.14,  3a/902  is  a  function  only  of  [ (tt- 02)  +  0  ]  as  expected, 

since  the  intensity  of  light  at  (r2,  @2) ,  as  given  by  Eq  7.26,  must 

be  independent  of  the  choice  of  coordinate  system.  We  therefore 
modify  Eq  7.26  to  read 

J(r2)/S  =  |  { Cot  +  01)/[r2(1r-02  +  0^]}  (3a/302)  |  7.28 
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Figure  16.  Description  of  the  stellar  disc  in  terms  of  angular  displace 
ments.  The  relationships  of  Eqs  7.23  and  7.24  are  easily 
obtained  by  noting  that  0Q  =  0cos<j>  +  0j  cos<t> j  nd  that 

0  s in<J>  =  4> j  sin4>j.  These  two  equations  can  be  solved  for 


0,  and  <J>. . 


We  are  now  in  a  position  to  discuss  what  an  observer  should  see  as 
two  stars  are  approximately  aligned.  Eqs  7.23  and  7.24  permit  us  to 
pick  an  arbitrary  point  on  the  circumference  of  the  more  distant  star 
and  to  compute  its  angular  displacement  from  the  nearer  star  by  means 
of  Eq  7.22.  Equation  7.28  permits  us  to  compute  the  intensity  of  light 

from  that  point.  The  closest  binaries  to  be  resolved  optically7  are 

7.7  x  10  radians  apart,  and  the  limit  of  resolution  of  a  one  meter 

telescope  operating  at  0.5u  is  approximately  5  x  10  7  radians.  The 
nearest  star  is  approximately  four  light  years  away  and  the  angle 

-8 

subtended  at  that  distance  by  the  solar  diameter  is  3.4  x  10~  radians. 
Therefore  stellar  image  distortion  is  likely  to  be  considerably  beyond 
our  observational  capability.  The  three  images  expected  from  a  nearly 
aligned  pair  of  stars  are  depicted  in  Figure  17,  where  the  geometry 
is  described  by  Figure  15  and  r^  and  ^  are  four  light  years  in  distance. 

The  two  stars  are  assumed  to  be  like  the  sun  in  size,  mass  and  light 
emission. 

The  vertical  positions  of  the  circles  show  the  apparent  angular 
separations  of  the  three  images  and  the  size  of  the  circles  show  the 
relative  brightness.  The  image  of  the  nearer  star  is  assumed  to  be  on 
the  horizontal  zero  line.  It  is  seen  that,  although  three  images  are 
theoretically  lossible  for  two  nearly  aligned  stars,  the  image  from 
the  more  distant  star,  corresponding  to  the  closer  approach  and  greater 
bending  of  the  light  rays,  quickly  becomes  very  dim  as  the  stars  move 
out  of  alignment.  It  also  moves  closer  to  the  image  of  the  nearer  star. 
The  other  image  from  the  more  distant  star  quickly  becomes  "normal"  in 
brightness  and  angular  displacement  as  the  stars  move  out  of  alignment. 

The  most  likely  observation  of  gravitational  focusing  is  likely 
to  be  associated  with  spectroscopic  binaries,  where  alignment  should 
be  coincident  with  the  minimum  Doppler  shift.  The  optimum  conditions 
for  observation  would  be  to  have  the  observer  in  the  plane  of  revolution 
of  two  stars  which  are  coupled  closely  enough  gravitationally  to 
have  a  reasonably  short  period  of  revolution.  One  star  should  be  big 
and  bright  and  the  other  should  be  more  massive  but  very  small  to 
avoid  occlusion  of  too  much  light.  The  requirements  for  alignment,  close 
association,  and  differences  in  type  of  star  will  of  course  eliminate 
most  observed  spectral  binaries,  but  it  would  be  interesting  to  search 
for  one  which  would  verify  gravitational  focusing. 


7 

E.  A.  Fath,  The  Elements  of  Astronomy,  4th  Edition,  McGraw-Hill  Book  Co., 
New  York  &  London,  1944,  p  284. 
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Massive  galaries  with  dense  star  clusters  at  the  center  should 
show  limb  darkening  and  center  brightening  because  of  gravitational 
light  focussing,  but  without  resolution  of  individual  stars,  the 
effect  could  not  be  attributed  with  certainty  to  gravitational 
focussing.  Likewise,  stars  approaching  a  black  hole  condition  should 
be  bright  in  the  center  and  dark  around  the  edges. 
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Procedure 


Choose  r,9,  and  U,  l<r<°°,  O<0<2w,  0<U_<°°. 

Compute  |a+|  and  |a  |.  (Eq  2). 

Compute  0'  (Eqs  5  and  7)  for  |a+|  and  |a  |.  These  values  are 

labeled  0'  and  d' 

+ 

Compute  (9a/39)  +  and  (3a/ 90)  (Eq  3).  Use  0  for  (3a/90)+  and  (2ir-0) 

place  of  0  in  calculating  (3a/ 30)  . 

Compute  (r-t)+  and  (r*t)  (Eq  10).  Use  0  for  (r-t)+  and  (2tt-0)  for 
in  (r*t)  . 

Compute  a^  (Eq  12). 

Compute  |j(r,0,a+)|  from  Eq  8  using  |a+|,  (3a/30)+  and  (r*t)+ 
as  calculated  above.  Set  |j(r,9,a+)|  equal  to  zero  for  0>0'+  and  a+ 

Compute  |j(r, 0,a_)|  from  Eq  8  using  [a  [,  (3a/30)  and  (r*t)  as 
calculated  above.  Set  |j)r, 0,a  )|  equal  to  zero 

Compute  <f>(r,  0)  from  Eq  1. 

If  r«t  =  0  above,  use  Eq  4  and  11  in  Eq  9,  using 
procedure  as  before. 


for  a  <a  . 
-  g 


the  same  general 


in 


0 


<a  . 
-  g 
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APPENDIX  2 


Computation  of  Apparent  Stellar  Brightness  and  Image 
Distortion  Caused  by  Gravitational  (Schwarzschild)  Bending  of  Light  Rays 


The  physical  constants  assumed  for  the  calculation  are  listed  in 
Table  A2-1.  From  Eq  7.1  and  the  constants  of  Table  A2-1  we  see  that 

-8 

a  grazing  orbit  is  bent  by  an  amount  of  9  x  10"  radians.  The 

-8 

subtended  angle  of  the  nearer  star  is  3.35  x  10  radians.  The 
formulae  used  for  the  calculation  are  listed  in  Table  A2-2.  The 
procedures  for  the  calculation  are  listed  in  Table  A2-3. 
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Table  A2-1 

Physical  Constants  Assumed  for  the  Calculation 


G  =  6.67  x  10'11  (S.I.) 

O 

C  =  3  x  10  m/sec 

M  =  2  x  1030  kg 

rl  =  r2  =  4  x  10  m 
k  =  4GM/C2  =  6000  m 
R  =  6.7  x  108  m 


Universal  Gravitational  Constant 
Velocity  of  light 

Mass  of  sun 

Distances  from  star  1  to 
star  2  to  observer 

Radius  of  Sun 
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Table  A2-2 


Formulae  Used  for  Computations 

a  +  =  fr/4)  +  l(n-62)2  +  8k/r]1/2}  2 

a  1  =  fa  1  )r2  2 

'  afj  =  fl/4)  ±  f1/-*)  [C^-e2)]/[(ir-e2)2  +  8k/r]1/2  3 

*♦  =  k/a.  4 

6+  =  +  K  ~  \  5 

J+  w  r2/s  =  |[o+/(ir-e2)]Oa+/ae2)|  6 
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Table  A2-3 

Procedures  for  Calculations 


1. 


2. 

3. 

4. 

5. 

6. 

7. 


Compute  a+  and  a-  from  Eq  1.,  0  <  tt- 0^  <  10 

Compute  a+  and  a-  from  Eq  2,  using  a+  to  calculate  a+,  etc. 

Compute  |C3a/362)+|  and  |(3a/3e2)|-  from  Eq  3. 

Compute  \p+  and  i|i-  from  Eq  4. 

Compute  6+  and  6-  from  Eq  5. 

2  2 

Compute  J+(r)r  /S  and  J-(r)r  /S  from  Eq  6. 

Plot  6+  and  6-  versus  ir-Q^  a^d  label  each  point  with  the  value 
2  2 

of  J+(r)r  /S  and  J-(r)r  /S  as  appropriate. 


64 


I 


DISTRIBUTION  LIST 


No.  of 

Copies  Organization 

12  Commander 

Defense  Documentation  Center 
ATTN:  DDC-DDA 
Cameron  Station 
Alexandria,  VA  22314 

1  Commander 

US  Army  Materiel  Development 
and  Readiness  Command 
ATTN:  DRCDMD-ST 
5001  Eisenhower  Avenue 
Alexandria,  VA  22333 

2  Commander 

US  Army  Armament  Research 
and  Development  Command 
ATTN:  DRDAR-TSS 
Dover,  NJ  07801 

1  Commander 

US  Army  Armament  Materiel 
Readiness  Command 
ATTN:  DRSAR-LEP-L,  Tech  Lib 
Rock  Island,  IL  61299 

1  Commander 

US  Army  Aviation  Research 
and  Development  Command 
ATTN:  DRSAV-E 
12th  and  Spruce  Streets 
St.  Louis,  MO  63166 

1  Director 

US  Army  Air  Mobility  Research 
and  Development  Laboratory 
Ames  Research  Center 
Moffett  Field,  CA  94035 

1  Commander 

US  Army  Communications  Rsch 
and  Development  Command 
ATTN:  DRDC0-PPA-SA 
Fort  Monmouth,  NJ  07703 


No.  of 

Copies  Organization 

1  Commander 

US  Army  Electronics  Research 
and  Development  Command 
Technical  Support  Activity 
ATTN:  DELSD-L 
Fort  Monmouth,  NJ  07703 

2  Commander 

US  Army  Missile  Research 
and  Development  Command 
ATTN:  DRDMI-R 

DRDMI-YDL 

Redstone  Arsenal,  AL  35809 
1  Commander 

US  Army  Tank  Automotive  Rsch 
and  Development  Command 
ATTN:  DRDTA-UL 
Warren,  MI  48090 

1  Director 

US  Army  TRADOC  Systems 
Analysis  Activity 
ATTN:  ATAA-SL,  Tech  Lib 
White  Sands  Missile  Range 
NM  88002 

Aberdeen  Proving  Ground 

Dir,  USAMSAA 
ATTN:  DRXSY 

DRXSY-MP ,  H.  Cohen 
Cdr ,  USATECOM 

ATTN:  DRSTE-TO-F 
Cdr,  USACSL,  EA 

ATTN:  ATSL-CLC,  Bldg  E1570 
Dr.  Tom  Welch 

Dir,  Wpns  Sys  Concepts  Team 
Bldg.  E3516 ,  EA 
ATTN :  DRDAR-ACW 


65 


USER  EVALUATION  OF  REPORT 


Please  take  a  few  minutes  to  answer  the  questions  below;  tear  out 
this  sheet  and  return  it  to  Director,  US  Army  Ballistic  Research 
Laboratory,  ARRADCOM,  ATTN:  DRDAR-TSB,  Aberdeen  Proving  Ground, 
Maryland  21005.  Your  comments  will  provide  us  with  information 
for  improving  future  reports. 

1 .  BRL  Report  Number^ _ 

2.  Does  this  report  satisfy  a  need?  (Comment  on  purpose,  related 
project,  or  other  area  of  interest  for  which  report  will  be  used.) 


3.  How,  specifically,  is  the  report  being  used?  (Information 
source,  design  data  or  procedure,  management  procedure,  source  of 
ideas,  etc.)  _ 


4.  Has  the  information  in  this  report  led  to  any  quantitative 
savings  as  far  as  man -hours/ contract  dollars  saved,  operating  costs 
avoided,  efficiencies  achieved,  etc.?  If  so,  please  elaborate. 


5.  General  Comments  (Indicate  what  you  think  should  be  changed  to 
make  this  report  and  future  reports  of  this  type  more  responsive 
to  your  needs,  more  usable,  improve  readability,  etc.) _ 


6.  If  you  would  like  to  be  contacted  by  the  personnel  who  prepared 
this  report  to  raise  specific  questions  or  discuss  the  topic, 
please  fill  in  the  following  information. 

Name : 


Telephone  Number: 
Organization  Address: 


